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Abstract. We investigate the accuracy of Eulerian perturbation theory for describing the matter 
and galaxy power spectra in real and redshift space in light of future observational probes for precision 
cosmology. Comparing the analytical results with a large suite of N-body simulations (160 independent 
boxes of 13.8 (Gpc/ft.) 3 volume each, which are publicly available), we find that re-summing terms in 
the standard perturbative approach predicts the real-space matter power spectrum with an accuracy 
of < 2% for k < 0.20/i/Mpc at redshifts z < 1.5. This is obtained following the widespread technique 
of writing the resummed propagator in terms of 1-loop contributions. We show that the accuracy of 
this scheme increases by considering higher-order terms in the resummed propagator. By combining 
resummed perturbation theories with several models for the mappings from real to redshift space 
discussed in the literature, the multipoles of the dark-matter power spectrum can be described with 
sub-percent deviations from N-body results for k < 0.15/i/Mpc at z < 1. As a consequence, the 
logarithmic growth rate, /, can be recovered with sub-percent accuracy on these scales. Extending 
the models to massive dark-matter haloes in redshift space, our results describe the monopole term 
from N-body data within 2% accuracy for scales k < 0.15/i/Mpc at z < 0.5; here / can be recovered 
within < 5% when the halo bias is known. We conclude that these techniques are suitable to extract 
cosmological information from future galaxy surveys. 
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1 Introduction 

Galaxy clustering is a key observational probe to study the large-scale structure of the Universe. 
The shape of the galaxy power spectrum, bispectrum and higher-order moments, contain information 
about the matter content of the Universe, about gravity and also about possible non-Gaussian initial 
conditions. In galaxy surveys, the galaxy distribution is distorted along the line of sight due to peculiar 
velocities that cause Doppler shifts, namely redshift space distortions (RSD). As these distortions 
depend on the growth of structure, they offer a complementary technique (to studies of the cosmic 
expansion history) to measure the matter content or to test gravity e.g., [1-3]. On very large scales 
and at higher redshifts, the RSD can be described by linear theory. However, on smaller scales and 
at later epochs, non-linearities start to play an important role and must be taken into account to 
accurately estimate the cosmological parameters from observational data. 

Current surveys like BOSS 1 [4], and future missions like EUCLID 2 [5] will soon provide un- 
precedented datasets about the distribution of galaxies on large scales. In order to extract useful 
information from data of this quality we need more accurate theoretical models of structure forma- 
tion. Standard perturbation theory (SPT) is the straightforward way of proceeding and has been the 
workhorse in the field for decades. However, its practical applications to statistics of the density field 
provide limited accuracy in both real and redshift space. A number of studies have investigated the 
fidelity of different theoretical models to describe the redshift space distortions and the possibility 
of extracting cosmological information from them [6-21]. This work explores the potential of de- 
scribing the redshift-space distortions combining different mappings between real and redshift space 
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with resummed perturbation theories. This is the basis to model the redshift-space distorted power 
spectrum of dark matter and dark-matter haloes. We then focus on the systematic errors that these 
models produce when they recover the logarithmic growth rate. For future surveys with forecasted 
statistical error on this quantity at the % level, the accuracy of the analytic description of clustering 
must be such that residual systematic errors due to modeling is well below this level. A fast, analytic 
description of clustering, calibrated on N-body simulations is an approach fully complementary to one 
based entirely on N-body simulations e.g. ,[22, 23]. 

In particular, we start by describing the real-space power spectrum using 1- and 2-loop standard 
perturbation theory and the so-called renormalized (or resummed) perturbation theory [24, 25]. We 
combine these models with different methods to obtain the redshift-space power spectrum: i) the 
Kaiser model [26], ii) the Scoccimarro model [27] and iii) the Taruya et al. model [28] (TNS model, 
from now on). For dark matter, we additionally account for the Finger-of-God (FoG) effect produced 
by virial motions on small scales by introducing a phenomenological damping term. This term is not 
needed in the case of haloes, as we only consider isolated haloes which are not part of a larger host halo. 
We compare these different methods with the results of a large suite of N-body simulations, focusing 
on the multipole expansion of the redshift-space power spectrum. Our suite of N-body simulations 
sums up to a volume of 2212 (Gpc//i) 3 , which is much larger than the volume surveyed by any current 
or planned experiment, ensuring that statistical errors in the simulations are negligible. 

This paper is organized as follows. In §2 we present the basic theory of redshift-space distortions. 
In particular, in §2.1 we review standard perturbation theory and renormalized perturbation theory, 
extending the latter to higher-order propagators than previously considered. In §2.2 we present 
different RSD models while Finger-of-God effects are discussed in §2.3. In §3 we provide details of the 
N-body simulations and describe the halo catalogues used in this paper. In §4 we compare the results 
of our models for the real-space power spectrum and compare them with N-body simulations for the 
dark-matter case. We also consider the RSD models mentioned above, focusing on their capacity to 
recover the logarithmic growth rate /, both for dark matter and massive haloes. Finally in §5, we 
discuss and summarize the obtained results. Appendices A and B contain details about standard and 
resummed perturbation theory providing a justification of the formulae presented in §2.1. Note that 
our discussion of resummed theories does not assume a field-theory background and is supposed to 
be accessible to all readers. 



2 Theory 

The matter-matter real-space power spectrum Pss (k) , the Fourier transform of the two-point corre- 
lation function, is the simplest statistic of interest one can extract from the dark matter overdensity 
field in Fourier space, <5(k), 

(S(k)S(k')) ee (2^) 3 <5 D (k + k')P 5S (k) , (2.1) 

where 6 D denotes the Dirac delta function and (...) the ensemble average. Under the assumption 
of an isotropic Universe, the power spectrum in real space does not depend on the direction of the 
k-vector. Since we only have one observable Universe, under the hypothesis of ergodicity the average 
(...) can be taken over all different directions for each k-vector. 

The mapping between the radial coordinate in real space and the radial coordinate in redshift 
space is given by the Hubble flow and the Doppler effect due to the peculiar velocities, v. Since only 
the radial distance is computed from the measured redshift, the two angular coordinates remain the 
same in both real and redshift space. In this paper, we adopt the distant observer approximation, i.e., 
we assume that all line-of-sights are virtually parallel to each other. If we identify this direction with 
the third axis of our coordinate system, the mapping from real-space coordinates x to redshift-space 
coordinates s reads 

s = x + #rf* 3 ' ^ 

ti (a)a 

where H(a) is the Hubble parameter at the scale factor a and x.3 denotes the unit vector of the third 
axis. Using the scaled velocity field u = —v/[H(a)af(a)} where / is the logarithmic derivative of the 
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linear growth factor with respect to the scale factor, / = d In D(a)/d In a, we can write this mapping 

as 

s = x-/it 3 (x)x 3 . (2.3) 
The 2-point correlation function in redshift space is then defined by 

S(r) = <*(s + r)<J(8)) (2.4) 

and its Fourier transform, 

P| 5 (fc,/i) = J d 3 r^(r)cx P Hk- r), (2.5) 

is the power spectrum in redshift space. Here, /i = k • k 3 is the cosine of the angle between the vector 
k and the line-of-sight. The redshift space power spectrum, can also be computed from Eq. 2.2. In 
this mapping the mass must be conserved, which implies [1 + 5^ (s)]d 3 s = [1 + 5(r)]d 3 r. Thus the 
transformation from S(r) to <5^(s) at linear order in 8 and v reads, 

*•>(-)= *(x)-^g« (2.6) 

where V3 is short for The two-point correlation function in Fourier space, (<5^ s )(k)(5^ s ^(k')) is 

then [28], 

= / Ae* r {e- tt ^ A ^ 0(x) + /V 3U3 (x)][5(x') + /V 3 u 3 (x')]>, (2.7) 

where r = x — x' and Ait z = u z (x) — u z (x'). Note that we have written u instead of v to make 
the dependence on / explicit. In Eq. 2.7 the enhancement and damping effect of the redshift space 
distortions on the power spectrum arc manifest. The enhancement due to RSD, also known as Kaiser 
effect, is produced by the the +/V3U3 terms in Eq. 2.2, that increase the overdensity <5(x). These 
terms, represent the coherent distortions by the peculiar velocities along the line-of-sight direction, 
and are controlled by the growth factor parameter /. On the other hand, the damping effect comes 
from the exponential factor in Eq. 2.7. This term is mainly due to the small scale velocity dispersion 
around the most clustered regions, and produces the suppression of power at small scales in the power 
spectrum. 

2.1 Perturbation theory in real space 

In order to describe the non-linear matter power spectrum in redshift space we first need a theory 
that is able to provide an accurate description of the power spectrum in real space. There are 
several models that attempt to do this task: the halo model (see [29] for a review), HALOFIT [30], 
cosmological standard perturbation theory [31-40] , and other perturbation theories approaches based 
on Lagrangian perturbation theory [41-45], time renormalization [46, 47], Eulerian renormalized (or 
resummed) perturbation theories [24, 25, 48-52] and closure theory [53] (see [54] for comparison of 
some of these theories). 

In this paper we focus on two approaches: standard perturbation and resummed perturbation 
theory, both in Eulerian space. 

Standard perturbation theory (SPT hereafter) consists of expanding the statistics of interest 
as a sum of infinite terms, where every term correspond to a n-loop correction. For the power spectrum 
the SPT prediction is written as (see appendix A for the explicit formulae of SPT terms), 

P SPT (fc) = p(°)(fc) + P«(/fc) + P {2 \k) + .... (2.8) 

The 0-loop term correction is just the linear power spectrum, pW(fe) = P lm (fc). The 1-loop term is 
expressed as a sum of 2 different subterms, 

P (1 \k) = 2P 13 (k) + P 22 (k), (2.9) 
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where the subscripts i and j refer to the perturbative order of the terms SQs.) used in eq. 2.1 to compute 
the power spectrum Py(fc). In this case, both P13 and P 22 requires a 2-dimcnsional integration (after 
exploiting rotational invariance). The 2-loop term is the sum of three different subterms, 



P^(k) = 2P 15 (fc) + 2P 24 (fc) + P 33 (k). 



(2.10) 



In this case, all these three terms, require a 5-dimensional integration (after exploiting rotational 
invariance). One can keep going to higher-order terms. However, the 3-loop correction term requires 
the computation of 8-dimensional integrals and the 4-loop correction term 11-dimensional integration. 
For practical and computational reasons, one does not usually go beyond the 2-loop correction terms. 
Thus, one has to truncate Eq. 2.8 series at some loop order. Truncating it at pW and at P 1 - 2 ' term 
is what is respectively called in this paper, 1L-SPT and 2L-SPT. 

Renormalized perturbation theory (RPT hereafter) attempts to reorganize the perturbative 
series expansion of SPT and resum some of the terms into a function that can be factorized out of 
the series. This function is usually called the resummed propagator and we refer to it as Af. In order 
to make the resummation possible, all the kernels of the P^ terms have to be expressed as a product 
of kernels that correspond to full-mode-coupling terms and full-propagator terms (see appendix B for 
details). The full-mode-coupling kernels are those kernels contained in P nn {k) terms that contain a 
coupling of the form, k qi — . . . — q„_i (see Eq. B.5). The full-propagator kernels are those contained 
in P^ terms of the form P\ n with no- mode coupling term (see Eq. A. 4 and A. 6 as examples). The 
resulting expression of resumming terms in this way is (see appendix B and Refs. [24, 25] for a full 
derivation) , 



P 



RPT-Wi 



(k,z) 



P lin (k 7 z) + P 22 (k,z) + Pi 3 L (k,z) + 



+ Pt~ 1)L (Kz) + ..]M % {k) 



(2.11) 



where the term P„ 



(n-l)L 



is the part of the P nn term that describes a full-mode coupling (see Eq. 



B.5). In spite of the resummation, Eq. 2.11 contains an infinite series as Eq. 2.8 and has to be 
truncated after a certain number of loops. However, some of the infinite terms of Eq. 2.8 have now 
been reorganized into the Ni function. We will refer as lL-RPT-A/i and 2L-RPT-A/i the truncation 
of Eq. 2.11 at 1- and 2-loop, respectively. Finally, the form of the function A/i depends on the way 
we approximate the kernels in the resummation process. In the case the kernels are approximated 
according to the Zel'dovich approximation (see Eq. B.l), they are expressed as a product of 0-loop 
propagators and the resulting function Ao is, 



Ao(fc) = cxp 



1 



with, 



in 

T 



dq 



k 2 at 



P hn (q)- 



(2.12) 



(2.13) 



When the kernels are approximated as a product of 1-loop propagator kernels (see Eq. B.17-B.19) 
the resulting A/"i function is, 

AMfc^cxpfP^fcVP 1111 ^)]. (2.14) 

This is the expression presented in [56] (sec Eq. B.31 for the correspondence). Expressing the kernels 
as a product of 2-loop propagators, (see Eq. B.32-B.37) yields the function N 2 , 



N 2 {k) = cosh 



; 2Pi 5 (fc) 
P lin (fc) 



Pl 3 (fc) / P' in (fc) 

P lin (fc)V 2P 15 (fc) 



sinh 



^PisCfc) 
P lin (fc) 



(2.15) 



Note that the angular part of the propagator terms, P13, P15, P17, 3 . . .is analytically integrable for 
any shape of the power spectrum. Thus, the 2-dimensional integration of P13 can be reduced to a 



3 Pi 7 is required for the 3-loop expansion of the resummed propagator, A/3 (see Eq. B.43). 
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1-dimcnsional integration (see Eq. B.31). In the same way, the 5- and 8-dimensional integrations 
of the terms Pig and Pyr, are reducible to 2- and 3-dimcnsional integration. However this is a hard 
task due to the symmetrized kernels, which are constructed as a sum of 5! = 120 and 7! = 5040 
terms respectively (see Eq. A. 14). Because of that, we stop at 2-loop. In this paper, the Af- 
function is computed considering Pig as a 5-dimensional integral. If more accuracy is needed, the 
3-loop resummed propagator written in Eq. B.43 can be used. These extensions in the resummed 
propagator could be easily incorporated in the current public codes for the RPT [50, 54, 56]. 

2.2 Perturbation theory in redshift space 

In order to describe the non-linear matter power spectrum in redshift space (Eq. 2.5) we need a model 
that, given the non-linear power spectrum in real space, is able to "map it" to the power spectrum 
in redshift space. There are several models that attempt to do this task. While in principle the same 
perturbation theory approach used for the dark matter could be employed to model also the velocity 
and the density velocity coupling yielding therefore a real-to-redshift space mapping, it has become 
clear in the literature that the redshift space clustering and in particular the redshift space power 
spectrum is not well described perturbatively. In fact, highly non-linear scales are superimposed to 
linear scales by the real-to-redshift space mapping, as realized in the seminal papers by [26, 59]. Some 
of this effect is even visually apparent in the galaxy distribution as the so-called "Fingers of God" 
effect. In this paper we consider physically motivated, but phenomenological models, we study the 
Kaiser model [26], the Scoccimarro model[27] and the TNS model [28]. All these models propose a 
functional form of P s (k) that depends on real-space statistics. 

The simplest model, is the Kaiser model proposed by Nick Kaiser 25 years ago [26], 



where b(k) is a possibly scale-dependent biasing function, which relates the observable tracers to the 
underlying dark matter distribution. This expression is obtained when Eq. 2.7 is treated linearly. 
Because of that, in principle Pss should refer to the matter-matter linear power spectrum. However, 
the prescription P]' 5 n — > Pgg has been demonstrated to work better. With this recipe, the Kaiser 
model is usually known as 'non-linear Kaiser'. In this paper we refer to Eq. 2.16 with a non-linear 
Pss simply as Kaiser model. 

The Scoccimarro model proposes that the redshift-space power spectrum is given by [27], 



where Pse and Peg are the velocity-matter and velocity-velocity power spectra in real space, respectively, 
and are defined by 



where 6>(k) = [-ik ■ v{k)]/[af(a)H(a)]. Note that Eq. 2.17 tends to 2.16 when P S g and P ee tend to 
Pss- This is the case for the linear regime in SPT. 

The TNS model [28] takes into account the cross interaction due to linear and non-linear 
processes. This produces two extra terms to Eq. 2.17, 



P s (k, n) = [b{k) 2 Pss{k) + 2b(k)ii 2 fP S9 (k) + f^P ee {k) + A(k, (m, b) + B(k, ^, b)] . (2.20) 



where the A and B terms, arise from the interaction between the enhancement terms and the damping 
terms in Eq. 2.7. Also [57] have provided a complementary explanation and a possible generalization 
of the A and B terms. However, in this paper we use the basic formalism presented by [28], 



P s (k,fl) = (b(k) + f 2 [l) 2 Pss(k), 



(2.16) 



P s (k 7f i) = [b(k) 2 P ss (k)+2b(k)^fP se (k) + f 2 ^P d e(k)] , 



(2.17) 



(d(k)0(k')) = (2nfS D (k + k')Pse(k), 
(0(k)0(k')) = (27T) 3 5 D (k + k')Pee(k), 



(2.18) 
(2.19) 




(2.21) 



(2.22) 
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where, 



and 



F( q ! = § \ b(q)P 5 e(q) + f^Pee(q) } - ( 2.23 i 



!2,T) :i ^,(k 1; , :: iL , ,(k ! . k,. k,) = ( H(k, ! { b(k 2 )S(k 2 ) + /^0(k 2 )| {b(k 3 )6(k 3 ) + /|p(k 3 ; 



(2.24) 



with ki23 = ki + k.2 + k3. Since the functions A and B require an integration over the whole range of 
momenta q, we cannot use RPT predictions to compute them. Furthermore, A requires the crossed 
bispectra between S and 9, whose computation at 1- and 2-loop requires more effort than in the power 
spectrum case. Since we expect A and B to be small compared to Pss, Pso and Pee [28], in this paper 
we compute A and B assuming the leading terms for the power spectrum and bispectrum inside the 
integrals of Eq. 2.21 and 2.22. For the rest of the power spectrum terms of Eq. 2.16, 2.17 and 2.20 
we use the perturbative approaches described in section 2.1. 

None of these models is able to account for non-linear effects such as Fingers of God. These 
effects have to be included ad hoc through a function that damps the power spectrum in redshift 
space at small scales. 

2.3 Fingers of God 

The effects of small-scale velocities are not completely included in the models presented in section 
2.2. Both the Kaiser and Scoccimarro model ignore them as we have already commented. Only the 
TNS model takes them into account, but only as a cross term with large-scale squashing. Thus, the 
effect of these small-scales velocities has to be inserted as a multiplicative damping function into the 
models described in section 2.2 [58-60], 

P s (k, M ; z) -+ P s (k, fx; z)Dl oG {k, W z). (2.25) 

The most used prescriptions for that, are the Gaussian and the Lorentzian functions that both depend 
on a rcdshift-dependent parameter, o~(z) = o~qD(z), 

Dfe(*,M, *;<*)= 1 + . 5[kf J oD{z)f{z)] 2 Lorentzian, (2.26) 
D^(k,fi,z;a ) = exp{-0.5[kno- D(z)f(z)} 2 } Gaussian. (2.27) 

where cr = a(z = 0). Theoretically it has been suggested [27] that <j(z) can be computed analytically 

as, 

2, s 4n f dq 



v\z) = — J j^Pee(q,z). (2.28) 

Such a parameter is physically motivated as it tries to enclose the effect of the velocity dispersion 
of dark matter particles. Similar modeling was discussed e.g. in [61] and in [27], although the two 
have different interpretations. In any case, the fit and theoretical numerical value for a need not to 
coincide with the actual value of the velocity dispersion. The theoretical value is computed in the 
linear approximation and the fit value is obtained using a Gaussian approximation while the actual 
velocity dispersion is highly non-gaussian even on large scales. More discussion on this in [27]. In 
most of the cases, this modeling has shown a very poor agreement with N-body simulation results. 
Therefore, in this paper we always treat uq as a free parameter to be fit from N-body simulations. 

3 Simulations 

The simulations used in this paper model the structure formation on very large scales within a flat 
ACDM cosmology consistent with current observational data. The adopted cosmological parameters 
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are: Q A = 0.73, O m = 0.27 h = 0.7, n b h 2 = 0.023, n s = 0.95 and a 8 {z = 0) 0.8. Our suite of 
simulations consists of 160 independent runs with 768 3 particles in a box of length 2.4 h Gpc. Hence, 
each box has a volume of Vibox = 13.8 (Gpc//i) 3 and in total, we simulate about 2, 200 (Gpc//i) 3 . 

The initial conditions were generated at the starting redshift z = 19 by displacing the particles 
according to the second-order Lagrangian PT from their initial grid points. The initial power spectrum 
of the density fluctuations was computed with CAMB [62] . The simulations were performed with the 
GADGET-2 code [63] taking only the gravitational interaction into account. 

In this paper, we consider snapshots at z = 0, 0.5, 1, 1.5. In order to obtain the dark-matter 
field from particles we discretize each box using 512 3 grid cells. Thus the size of the Cartesian mesh 
is 4.68 h~ 1 Mpc. We assign mass to the cells using the cloud-in-cell prescription. Using a much 
higher resolution simulation, we checked that the power spectrum derived from the simulation data 
is accurate at the 1% level up to k < 0.2 /iMpc -1 . 

To identify the dark-matter haloes, we used the Amiga Halo Finder [64, 65], which defines a 
halo by the bound dark-matter particles inside a spherical overdensity equal to the so-called virial 
overdensity. We only consider haloes which are at least resolved by 40 particle. This leads to a 
minimum halo mass of 10 14 M Q //i. 

The errors associated to the statistical quantities measured from the simulations are obtained 
from the dispersion among the 160 runs: we report the error on the mean of the independent runs. 
We make the dark matter and halo power spectra publicly available and also provide the multipoles 
measured from the simulations used in this paper for possible comparisons 4 . 

4 Results 

In this section we compare and test the theoretical formalism presented in §2 (for both real and 
redshift space) with the N-body simulation data. We start with perturbation theory in real space, 
and we follow the theoretical predictions of different models for redshift-space power spectrum: Kaiser 
(Eq. 2.16), Scoccimarro (Eq. 2.17) and TNS model (Eq. 2.20). In particular, we consider here 1- and 
2-loop SPT, and 2-loop RPT with Afi and A/2- We focus on the multipole prediction of these models 
but also on the accuracy on recovering the logarithmic growth rate /. We explore this for both dark 
matter and for massive dark matter haloes. 

4.1 Performance of the perturbation theory approach: comparison to N-body simula- 
tions 

In this subsection we test the formalism presented in section §2.1, both SPT and RPT, and we compare 
them with the outcome of the N-body simulations. 

In Fig. 1-2 (left panels) we compare the power spectra obtained with different flavors of 
perturbation theory against the N-body simulation results at z = 0,0.5,1.0, and 1.5. Each power 
spectrum has been normalized to the linear non-wiggle power spectrum to reduce the dynamical 
range. In the left top subpanels the power spectrum of the N-body simulations is shown as black 
circles linked with a black dot-dashed line, the linear theory prediction is shown as a dotted black 
line, SPT in blue lines, RPT-A/o model in red lines, RPT-A/i model in green lines and RPT-A2 model 
in orange lines; where the infinite series of Eq. 2.8 and 2.11 are truncated at 1-loop (solid lines) and 
2-loop (dashed lines). In the bottom left subpanel the ratio of these models with N-body simulation 
data is shown with the same color and line notation. In that case, circle symbols correspond to linear 
prediction, square symbols to 1-loop truncation and triangle symbols to 2-loop truncation. 

We see that SPT with 1-loop truncation over-predicts the N-body power spectrum at all redshifts, 
and can only make an accurate prediction of the power spectrum (< 1% deviation) at very large 
scales: k < 0.05/i/Mpc for z — and k < 0.10/i/Mpc for z = 1. Going to 2-loop correction improves 
considerably the behavior of SPT, but one has to deal with the terms P33, P24 and Pi 5 that require 
a 5-dimensional integration. In this case, 2-loop SPT makes a very good prediction of the N-body 

4 http : //ice .ub . edu/~hector/Hector_Gil_Marin/Public .html 
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Figure 1. Left top subpanel: power spectrum at z = normalized to the linear non- wiggle power spectrum 
to reduce the dynamic range. N-body data in black circles and dot-dashed line. Different theoretical models 
are also shown: linear prediction (black dotted lines), SPT model (blue lines), RPT-Ao model (red lines), 
RPT-A/i model (green lines) and RPT-A/2 model (orange lines); for 1-loop truncation (solid lines) and 2-loop 
truncation (dashed lines). Left bottom subpanel: ratio between the power spectrum of N-body simulation 
and different PT models with the same color notation. Right top subpanel: Ao(fc) 2 (red solid line), A/i(fe) 2 
(green dashed line) and Afa(k) 2 (orange dotted line) at z = 0. Right bottom subpanel: <j 2 (red solid line), 
-2Pi 3 (fc)/[fc 2 P lin (fc)] (green dashed line) and -2 ln[A/" 2 (fc)]/fc 2 (orange dotted line) in units of (Mpc//i) 2 . 

power spectrum at z = up to relatively small scales (k < 0.20 /i/Mpc), but at z = 0.5 starts over- 
predicting it, and at z = 1 and z = 1.5 the over-prediction is a few percent at intermediate scales 
(* > 0.1/i/Mpc). 

The RPT-Ao model for both 1- and 2-loop truncation behaves accurately only at very large 
scales where it shows < 1% deviation respect to N-body simulations both at all redshifts, but breaks 
down at relatively large scales: k ~ 0.03/i/Mpc for 1-loop at z = and k ~ 0.07/i/Mpc for 2-loop at 

3 = 0. 

The RPT-A/i model presents at large scales a very good agreement with the N-body simulation 
results with the advantage that the breakdown happens at smaller scales: k ~ 0.1/i/Mpc for 1- 
loop and k ~ 0.15 for 2-loop at z — 0; at k ~ 0.15/i/Mpc for 1-loop and at k ~ 0.25 for 2-loop at 
z — 1. However, the for 2-loop truncation over-predicts the N-body power spectrum with a systematic 
~ 2% deviation at intermediate scales at all redshifts. This is due to the limitation in the resummcd 
propagator A/i, which is expressed as 1-loop expansion terms (in terms of P13). 

Finally, the RPT-A/2 model presents a modest improvement over previous models on large scales 
at z = 0, but breaks down already at k ~ 0.07/i/Mpc and k ~ 0.10/i/Mpc for 1-loop and 2-loop, 
respectively. However at z > 0, this model presents a better behavior, with < 2% accuracy, breaking 
down at k ~ 0.20/i/Mpc at z = 1 and fixing the systematic ~ 2% over-prediction observed for the 
RPT-A/i model. The accuracy of this model at z > 0.5 is then better than the 2-loop RPT-A/i and 
2-loop SPT. 

In general, we see that the model that best describes the N-body data at z < 1.5 is the RPT-A/2 
model at 2-loop truncation. The 2L-RPT-Ai model shows also good results but presents a ~ 2%- 
systematic over-prediction at intermediate scales. However at z = 0, RPT-A/i is able to reach smaller 
scales than RPT-A2, which breaks down at relatively large scales. Indeed, 2L-RPT-A2 works better 
than 2L-RPT-A"i and 2L-SPT at z > 0.5 but not at z = 0. 

In Fig. 1-2 (right panels) the behavior of the damping functions of RPT-A/i models is shown for 
the same redshift range. In particular, in the right top subpanel we show the scale dependence of the 
damping functions of Eq. 2.11: Af (k) 2 = enp(~k 2 a 2 ) (red solid line), Ni{k) 2 = exp [2Pi 3 (fc)/P lin (£;)] 
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Figure 2. Same notation that in Fig. 1 but for z = 0.5 (top panels), z = 1.0 (middle panels) and z = 1.5 
(bottom panels). 
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(green dashed line) and A/^fc) 2 ( see Eq. 2.15) (orange dotted line). We see that at large scales, all 
Mi functions converge to 1, and it is at intermediate and smaller scales (k > 0.03 h/Mpc) where 
the differences between these three models are significant. For the scales of interest, we see that 
Mo(k) < M%(k) < M\(k). Thus, as we include higher-order propagators in the computation of Mi, 
these functions oscillate about the 'true' damping function. This approach explains why RPT-Ao 
under-predicts N-body data and why RPT-A/i slightly over-predicts it. Every loop correction in the 
resummed propagator tends to the true value, but in a oscillatory way. For z > 0.5, the 2-loop 
correction in the resummed propagator seems to be sufficient for a < 1% prediction. However at 
z = 0, the convergence in the resummed propagator is still not reached for the 2-loop correction in 
M. In that case, higher-order loop corrections would be necessary to reach the < 1% deviation. In 
the right bottom panel we show the effective <r 2 , i.e. — 21n[A/i(fc)]/fc 2 , for the different orders in the 
resummed propagator: 0-loop (red solid line), 1-loop (blue dashed line) and 2-loop (orange dotted 
line). In that case, we see that at large scales the three models diverge, whereas at small scales all 
models seem to converge (at least for z > 0.5). 

In other words, it may happen that a lower-order approximation appears to work better than a 
higher-order approximation which, in principle, should be more accurate. This is due to a fortuitous 
cancellation of the truncation errors in the propagator and in the damping function. This cancellation 
does not hold for all redshifts and/or all cosmologies. The performance of an analytical approximation 
scheme must be quantified looking at different redshifts (or different cosmologies). 

The formalism presented here deals with the F n kernels, that correspond to the J-field. However, 
this formalism is also perfectly valid for the computation of Pse and Pee , only changing appropriately 
the F n kernels by the G n kernels as in SPT. 

4.2 Dark matter multipoles 

In order to obtain information about the growth rate / from the redshift-space distortions, it is 
convenient to work with the expansion in Legendre moments, Pg, defined as, 

P t (k) = {21+1) f d^P s {k,n)L t (fx), (4.1) 
Jo 

where Li are the Legendre polynomials of order £. For the first three non- vanishing Pg, 

L (x) = 1, (4.2) 

L 2 (x) = i(3x 2 -l), (4.3) 

L 4 (x) = i(35x 4 - 30x 2 + 3). (4.4) 
8 

According to linear theory (Kaiser model with Pgg 1 ) only the monopole (£ = 0), the quadrupole (I = 2) 
and the hexadecapole (£ — 4) are different from 0. In that case, for an unbiased tracer these three 
moments read, 

P (k) = P^(fc)(l + ^/ + i/), (4.5) 

P 2 (k) =P Un (A)(|/ + ^/ a ), (4.6) 

P 4 (fc) = P lin (fc) (A/^ (4 7) 

Hence, knowing the dark matter power spectrum in both real and redshift space, one can directly 
measure the growth rate / from any of these multipoles. It is interesting to note that the ratio 
between any of these multipoles does not depend (at large scales) on the real-space power spectrum 
and the ratio tends to a constant that only depends on / when k — > 0. However, non-linearities 
produce deviations from these formulae. Depending on the ability of modeling the non-linearity in 
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the redshift space distortions, we will be able to use information from non-linear scales to estimate / 
with accuracy. 

In this section we focus on checking the quality of the different theoretical models in predicting 
the multipole power spectrum of dark matter in redshift space. We focus on the models described 
in section 2.2, using as real-space inputs, the PT-theory approaches described in section 2.1. Here, 
we assume that / is known and we only fit the FoG parameter ctq assuming a Lorentzian damping 
function (Eq. 2.26), although no significant difference is observed when a Gaussian damping function 
is assumed. We allow <j$ to depend on z and on /c max 5 and find the best-fit value by minimizing 



Here, the subscript "sims" refers to the simulations and "theo" to the theoretical models described 
above. a\ s ; ms (fc) is the variance of the multipoles computed from the 160 simulations and k is the 
minimum k considered, which is set by the size of the simulation box but has little effect on the final 
results. Note that we neglect for simplicity any covariance between different fc-bins, which is a good 
approximation for small k and broad bins. 

In Fig. 3, we show the measurements of the monopole (top panels), quadrupole (middle panels) 
and hexadecapole (bottom panels) from N-body simulations (black empty circles) in top sub-panels. 
In all cases, the multipoles have been normalized to the non-wiggle linear prediction to reduce the 
dynamical range. In the top subpanels, different theoretical predictions are shown: the Kaiser model 
(dotted lines), the Scoccimarro model (dashed lines) and the TNS model (solid lines). The chosen 
real-space power spectrum for each of these models is: linear prediction (black dotted lines), 1L-SPT 
(red lines), 2L-SPT (blue lines), 2L-RPT-A/i (green lines) and 2L-RPT-7V 2 (orange lines). 

In all cases, the TNS model combined with both SPT and RPT predictions is the model that 
describes best the N-body results. Using the TNS model it is possible to achieve < 1% accuracy for 
the monopole up to k < 0.12/i/Mpc at z = and k < 0.17/i/Mpc at z = 1. The models are also very 
accurate for the quadrupole: at z — we can describe N-body data up to scales of k = 0.12/i/Mpc 
and at z = 1 up to k = 0.30/i/Mpc with a deviation < 2%. For the hexadecapole, the agreement is 
more modest: at z — we can only achieve ~ 10% accuracy up to scales of k = 0.20 h/Mpc at z = 
and ~ 5% at z = 1. 

Both the Scoccimarro and Kaiser model provide a reasonably good approximation on large 
scales but both fail to give an accurate description on mildly non-linear scales where baryon acoustic 
oscillations (BAO) are located. The difference between the TNS model and the other models is more 
evident for the quadrupole and hexadecapole possibly suggesting that non-linearities become more 
important for higher-order multipoles. 

The imprint of the BAO in the multipoles is clearly visible: note that the Scoccimarro and Kaiser 
models slightly over-predict the BAO amplitude, especially for the quadrupole, while the TNS model 
does better, although a trend towards the under-prediction is observed. All models correctly predict 
the BAO location. These considerations might be relevant for recovering in an unbiased way the 
angular and radial BAO information (separately) from forthcoming surveys. 

We conclude that the TNS model with the RPT and SPT models studied here, has the ability 
of describing the redshift space power spectrum monopole and quadrupole at z = and z = 1 within 
1 — 2% for k < 0.2 and the hexadecapole within about ~ 5%. 

We do not observe a crucial difference between 1- and 2-loop SPT. Also, no significant difference 
between using A/i and A/a for 2-loop RPT is detected. This indicates that on these mildly non-linear 
scales at redshifts < 1.5 the accuracy of the modeling of redshift space distortions is more important 
than that of the non-linear evolution of the real-space dark matter power spectrum. 

Because of that, for simplicity we focus on 1-loop SPT and 2-loop RPT with A/i when measuring 
/ from dark-matter multipoles in the next subsection. 

In Fig. 4 we show the best-fit value for ctq corresponding to the fits shown in Fig. 3 using the 
same color notation for the different models. Additionally we show the theoretical predictions of Eq. 

5 fc max is the maximum k used for the fit 




(4.8) 
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Figure 3. Multipoles corresponding to dark matter power spectrum: monopole (top panels), quadrupole 
(middle panels) and hexadecapole (bottom panels), for z — (left panels) and z = 1 (right panels), where / is 
fixed to the true value, and oo is the only free parameter, fit to N-body data. In the upper subpanels the value 
of the corresponding multipole is shown, normalized to the linear, no-wiggles value to reduce the dynamical 
range. In the lower subpanels the ratio between the N-body simulation data and different perturbation theory 
predictions is shown. Dotted lines correspond to the Kaiser model, dashed lines to Scoccimarro model and 
solid lines to the TNS model. Different perturbation theory models are shown: linear prediction (black dotted 
lines), 1L-SPT (red lines), 2L-SPT (blue lines), 2L-RPT-A/i (green lines) and 2L-RPT-A/" 2 (orange lines). In 
bottom subpanels the 2% deviation is marked with black dot-dashed horizontal lines. 
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Figure 4. Best-fit values of ao as a function of fc max corresponding to multipoles shown in Fig. 3 with 
the same color notation. As indicated, top, middle and bottom panels stands for monopole, quadrupole and 
hexadecapole fits alone. Left and Right panels show the result at z = and z = 1 respectively. Also theoretical 
predictions for ao are shown according Eq. 2.28: with Pee as input using 1L-SPT prediction in solid black 
line and using Pn n as input in dot-dashed black line. Error-bars are not shown for the sake of clarity, but are 
negligible for k > 0.05 /i/Mpc. 



2.28 using as an input the 1L-SPT prediction for Pgg (solid black line) and Pi; n (dot-dashed line). As 
indicated, top, middle and bottom panels correspond to the monopole, quadrupole and hexadecapole, 
whereas left panels show the result for z — and right panels for z = 1. 

We note that at z = 1 all the models produce a best-fit ao which is close the the theoretical 
predictions, although an overestimate is observed for the monopole and underestimation for the hex- 
adecapole, being the quadrupole the case which is closer to the theoretical prediction. For z = 0, the 
discrepancy between theory and best-fit value is larger. We will analyze again the agreement between 
theory and best-fit <jq in the next section, when will allow / also to vary. 

4.3 Estimating / from dark matter multipoles 

In the last section we have shown that the TNS model was able to describe well the multipoles when 
one free parameter was allowed to vary in order to account for the FoG effect. In this section we 
want to check the ability of these models to recover / from the dark-matter field. In this case, we will 
allow both / and <7o to freely vary. As before, in order to find the / and <To best fit parameters, we 
minimize the \ 2 value. 

In Fig. 5 we show the obtained values for / (top subpanels) and gq (bottom subpanel) as a 
function of the maximum scale used in the fitting, namely k max . We show the results at different 
rcdshift: z = (left panels) and z = 1 (right panels). Top, middle and bottom panels show the 
derived values for / and Co for each of the multipoles: monopole, quadrupole and hexadecapole as 
indicated. As in Fig. 3, dotted lines stands for Kaiser model, dashed lines for Scoccimarro model 
and solid lines for TNS model. For simplicity, we only show the results corresponding to 1L-SPT (red 
lines) and 2L-RPT-7Vi (green lines). 2L-SPT and 2L-RPT-A/2 yield similar results. In top subpanels, 
horizontal solid black line shows the true value for /, whereas black dot-dashed horizontal lines show 
the 1% and 2% deviation, as labeled. In the bottom subpanels, the horizontal lines stands for the 
theoretical predictions of Co according to Eq. 2.28 when Pgg (blue line) and P lm (orange line) are used 
as inputs. In the case of Pgg we use 1L-SPT prediction. All the error bars correspond to 1 — a errors for 
/ and ao, and have been computed from the contour in the f-ao space that corresponds to Ax 2 = 2.3. 
Since our N-body sample consists of 160 realizations (~ 13 3 [Gpc/ft.] 3 volume), the dispersion in the 
measured monopole, quadrupole and hexadecapole is small. Therefore the corresponding error bars 
for the recovered parameters are also small, especially in the case of the monopole and quadrupole. 
However, one should be aware that these statistical errors in / and Uq are not comparable to the 
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Figure 5. Estimates for / (top subpanels) and <to (bottom subpanels) from the dark matter multipoles of 
the N-body data: monopole (top panels), quadrupole (middle panels) and hexadecapole (bottom panels); for 
2 = (left panels) and 2 = 1 (right panels). Results from different theoretical models are shown: dotted lines 
are Kaiser model, dashed lines Scoccimarro model and solid lines TNS model. Green lines are 2-loop-RPT-A/i 
and red lines 1-loop SPT. In top subpanels, the true value of / is represented in a horizontal black-solid line, 
whereas 1% and 2% deviations are shown in the horizontal black dot-dashed lines, as indicated. In the bottom 
subpanels the theoretical predictions for ao (Eq. 2.28) are shown in blue (with Pee as input using 1L-SPT 
prediction) and in orange (with Pn n as input). 
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expected errors from future surveys. These errors just provide information about the uncertainties 
and shortcomings of the models. Note that for most of the studied models, these errors are much 
smaller than the expected ones from any galaxy survey. 

From Fig. 5, we see that the TNS model using both 1L-SPT and 2L-RPT-A/"i is the only model 
able to recover the value of / to 1% accuracy: for the monopole up to k ~ 0.15/i/Mpc for z = 
and up to k ~ 0.20/i/Mpc for z = 1; for the quadrupole up to k ~ 0.20ft./Mpc for z = and up 
to k ~ 0.25/i/Mpc for z = 1; and for the hexadecapole up to k ~ 0.15ft./Mpc for z — and up to 
k ~ 0.20/i/Mpc for z = 1. However, in the case of the hexadecapole, the statistical errors are too 
large to be able to claim 1% accuracy of the model predictions. Also note the difference between 
the theoretical value of <j$ and the best-fit value obtained from the N-body data: at z = 1 for both 
monopole and quadrupole, the best-fit <7o value for TNS model + 2L-RPT-A/"i yields a very similar- 
result as Eq. 2.28. However at z = there is a large discrepancy between these two values. 

4.4 Halo biasing and stochasticity 

So far, we have been able to recover the / parameter with high accuracy using the dark matter 
multipoles. However, galaxy redshift surveys consist of galaxies residing in dark matter haloes, which 
are biased and stochastic tracers of the underlying dark matter field. Furthermore, since we only 
consider massive isolated haloes, we do not expect any FoG effects in the halo redshift-space statistics. 

In this paper, we model the biasing, i.e. the relation between the halo overdensity and the the 
dark matter density contrast as, 

5 h (k) = b(k)8(k)+e(k), (4.9) 

where b(k) is the scale-dependent bias function and e describes a stochastic field. The stochastic field 
e stands for any physical or statistical process that produces a non-deterministic relation between the 
dark matter and the halo field. This includes the shot noise due to the discrete nature of haloes 6 . For 
a Poisson process the shot noise is inversely proportional to the mean number density, namely nh, 

^Poisson = — ■ (4.10) 

However, the formalism used here allows for other stochastic processes. Assuming the e field to be 
uncorrelated with 5, namely, (Se) = 0, the bias function can be written as, 

(6 h (k)S(k'))S D (k + k')(27rr gnhjfc) , , 

° [ > (d(k)J(k'))<5 D (k + k')(27r)3 " P mm (fc)' ^ > 

where the second equality stands for the numerator and denominator independently. Here the sub- 
scripts "h" and "m" stand for haloes and dark matter respectively. The power spectrum of the e field 
can be computed combining Eq. 4.9 and 4.11, 

P ee (k) = (e(k)e(k'))5 D (k + k')(2nf = (S h (k)S h (k'))5 D (k + k')(2^) 3 - (4.12) 

Finally, we define the (shot noise free) halo-halo power spectrum Phh as, 

P hh (fc) = ( ( 5 h (k)4(k'))<5 D (k + k')(2 7 r) 3 - P ee (k). (4.13) 

Note that by this definition the following equality holds, 

P hh (k) = b(k)P mh (k) = b(k) 2 P mm (k). (4.14) 

Hence, in this biasing scheme, the bias functions that relate the halo-matter and halo-halo power 
spectra to the matter-matter power spectrum are the same. This is a reasonable approximation at 
least on large scales [66], where the bias becomes linear. 

In Fig. 6 we show the scale dependence of the halo bias (left panel) and the e-field power spectrum 
(right panel) measured from the halo catalogues at different redshifts. Remember that the minimum 
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Figure 6. Left panel: Halo bias for different redshifts snapshots with M cu t = 10 14 Mq//i. Right panel: e-field 
power spectrum normalized to the Poisson prediction, associated to haloes with M cut = 10 14 Mq /h. In both 
panels, red lines are for z = 0, blue lines are for z = 0.5, green lines are for z = 1 and orange lines are for 
Z = 1.5. 
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Table 1. Number density of haloes at different redshifts with a minimum mass of M cut = W 14 M Q /h. 
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Figure 7. Signal-to-noise for haloes with mass above M cut = lO 14 M //i for z = (red line), z = 0.5 (blue 
line), z = 1.0 (green line) and z = 1.5 (orange line). For reference, the values corresponding to 1, 0.5 and 0.1 
of the signal-to-noise value have been plotted in horizontal black dotted lines. 

halo mass is M cut = 1O 14 M //i, which ensures that all the haloes have at least ~ 40 particles. The 
number density of haloes with this mass cut at different redshifts is shown in Table 1. Fig. 6 shows 
that the bias increases with z and with k. The P ee increases with z and slightly decreases with k. For 
Mcut — 10 14 Mq//i, the P ee is sub-Poissonian at low redshifts, 65% at z — and 90% at = 0.5, but 
turns out to be very close to the Poissonian prediction at high redshifts, z = 1 and z = 1.5. 

6 In N-body simulations, dark matter particles represent also a discrete field. However, the number density of particles 
is large enough to render this effect negligible. 
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In Fig. 7 we show the signal-to-noise ratio (or Phh/Pee) for the haloes studied here at different 
rcdshifts: z = (red line), z — 0.5 (blue line), z = 1.0 (green line) and z — 1.5 (orange line). 
Horizontal dotted lines mark the values Phh/Pee — 1-0, 0.5 and 0.1 as a reference. In Table 2 the scale 
at which the signal-to-noise ratio reaches these values is written for the same z-snapshots. We see that 
only for z = and z = 0.5 the signal-to-noise ratio is above 0.5 at scales k < 0.1/iMpc, whereas for 
z = 1.0 is above 0.5 only at very large scales, k < 0.03 h/Mpc and never for z = 1.5. Conservatively, in 
this paper we consider only the scales where Phh/Pmm ^ 0.5 to be suitable for extracting information. 
Hence, we do not study the halo power spectra at z = 1.5 and z = 1.0 because only at very large 
scales (where the behavior is linear) the signal-to-noise ratio satisfies this condition. 



z 





0.5 


1.0 


1.5 


Phh/Pee 
Phh/Pee 
Phh/Pee 


= 1.0 
= 0.5 

= 0.1 


0.168 h/M-pc 


0.065 h/Mpc 
0.132 h/Mpc 


0.033 h/Mpc 
0.225 h/Mpc 


0.031 h/Mpc 



Table 2. Scale where the noise starts to be comparable to the signal: Phh/Pee = 1-0, 0.5, and 0.1 for haloes 
with M cnt = 10 14 M Q //i at z = 0, 0.5, 1.0 and 1.5. 



4.5 High-bias halo power spectrum multipoles 

The multipoles for the halo power spectrum are defined by Eq. 4.1 just changing P 1 ^ lm (fc,/z) by 
P^ h (k,/i). In this case, we define the halo-power spectrum in redshift space by 

P{ h (k,fi) = (^(k)^(k'))^(k + k')(2vr) 3 -P ££ (fc), (4.15) 

where we assume that P ee does not depend on /i. For the monopole term, the shot noise subtraction 
is important. However, since we are assuming that P ee does not depend on /i, it is irrelevant for 
higher-order moments. Indeed, a /i-independent offset on P^ h (k,fj.) has no effect in the quadrupole 
and hexadecapole, only in the monopole. As in Eq. 4.5-4.7, at very large scales, the halo multipoles 
are written as, 

Po(k) = P lin (fc) (V) 2 + h(k)f + i/ 2 ) , (4.16) 
P 2 (k) = P lin (fe)(^W/+~/ 2 ), (4.17) 
P 4 (fc) = P lin (fc) (~/ 2 ) . (4.18) 

Note that for the case of biased tracer, / and b are degenerate in the linear regime when we treat 
b 2 P hn as input or when any ratio between different multipoles is used to constrain /. In this case, 
only the ratio among them, (3(k) = f /b(k) can be measured. However, as for the dark matter case, 
non-linearities cause deviations from these formulae. Depending on the ability of modeling the redshift 
space distortions, we will be able to use information from non-linear scales to estimate / with accuracy 
but also, in the case of biased tracers, we might be able to break the degeneracy between the bias 
and /. For both the Kaiser and the Scoccimarro model, / and b always appear in the /3 combination, 
only for the TNS model this degeneracy is not exact. 

Fig. 8 shows the halo monopole (top panels) and halo quadrupole (bottom panels) for z — 
(left panels) and z = 0.5 (right panels). Both / and b(k) have been set to their true values. The color 
and line notation is the same as in Fig. 3: dashed lines stands for Scoccimarro model and solid lines 
for TNS model. The different real-space power spectrum inputs are: linear prediction (black dotted 
lines), 1L-SPT (red lines), 2L-SPT (blue lines), 2L-RPT-A/"i (green lines) and 2L-RPT-A/" 2 (orange 
lines). Since we do not expect Fingers of God for isolated haloes, the theoretical models have no free 
parameters in this case. 
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Figure 8. Multipoles corresponding the halo-halo power spectrum: monopole (top panels), quadrupole 
(bottom panels), for z — (left panels) and z = 0.5 (right panels). In the upper subpanels the value of 
the corresponding multipole is shown, normalized to the corresponding non-wiggle linear prediction to reduce 
the dynamical range. In the lower subpanels the ratio between the N-body simulation data and different 
PT-predictions is shown. Dashed lines corresponds to Scoccimarro model and solid lines to the TNS model. 
Different PT models are shown: linear prediction (black dotted lines), 1L-SPT (red lines), 2L-SPT (blue lines), 
2L-RPT-A/1 (green lines) and 2L-RPT-A/2 (orange lines). In bottom subpanels the 2% deviation is marked 
with black dot-dashed horizontal lines. Vertical dot-dashed lines mark the regions where Phh/Pee — 1-0 and 
Phh/Pee = 0.5 as labeled. 



From Fig. 8, we see considerably different results compared to the dark matter case (see Fig. 
3). The accuracy of the modeling is reduced, especially for the quadrupole at z = 0. In the case of 
the halo-halo monopole, we see that the different PT models make very different predictions, while 
this was not the case for the dark matter. In particular we see that TNS + 2L-RPT-A2 is the only 
model able to make sub-percent predictions at z = up to k — 0.15/i/Mpc and at z = 0.5 up to 
k = 0.10/j/Mpc for the monopole. Any other PT theory + TNS yields worse results. We also see 
that Scoccimarro + 1L-SPT and Scoccimarro + 2L-RPT-A/i provide a good description at z — but 
not at z = 0.5. This is due to an (accidental) cancellation of two terms that go in opposite directions. 
The Scoccimarro model does not take into account the A and B functions of the TNS model, that 
add a positive contribution to P s (see Fig. 11 in the next section). On the other hand, in Fig. 1 we 
have seen that 1L-SPT and 2L-RPT-A/1 over-predict the true value for Pss at z = (but also for P$g 
and Pee)- For the dark matter power spectrum, both effects are small and negligible. However, for 
massive haloes, the high bias increases these effects. At z = 0.5, the TNS approach is clearly modeling 
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better the monopole than the Scoccimarro model, which under-predicts the N-body data by ~ 10%. 

In the case of the halo-halo quadrupole, we see that all the PT models make similar predictions, 
and that the main difference is due to the RSD model chosen but the modeling breaks down at rela- 
tively large, almost linear, scales. We see that the TNS model is able to describe well the quadrupole 
only up to k = 0.05 ft,/Mpc at z = and up to k = 0.10 h/Mpc at z = 0.5. 

In general we see that for the monopole, all the models describe the N-body data better at z = 
than at z = 0.5. This seems counter-intuitive, because at higher redshifts, non-linearities are less 
important and perturbation theory should work better. However, it could be explained by the fact 
that at z = 0.5 the signal-to-noise ratio is considerably less than at z = (see Table 2 and Fig. 7). 
For instance, recall that at z = 0, Phh/Pee = 1-0 at k = 0.168 h/Mpc whereas at z = 0.5 this happens 
already at k = 0.065 h/Mpc. Also the bias of the selected haloes grows with rcdshift, making more 
apparent any systematic error in the model (of both RSD and biasing). It is reasonable to expect 
(but it remains to be tested) that the model performance improves for lower mass -thus less rare and 
less biased- haloes. 

4.6 Simultaneously estimating / and b from halo multipoles 

In this section we show how well the / parameter can be recovered from the halo-monopole N-body 
data. Since in the last section we have seen that none of the models studied here is able to reproduce 
the quadrupole data for haloes with sufficient accuracy at the mildly non-linear scales we are interested 
in, we do not try to recover / from Pi. Instead we focus on the degeneracy between / and the bias 
in the monopole. According to Eq. 4.5, these two parameters are perfectly degenerate when P(k)b 2 
is set from observations or when the P2/P0 ratio is used to compute / and b. This is the case for the 
Scoccimarro and Kaiser models without Finger-of-God effects. However, non-linear terms, namely 
the A and B functions of the TNS model, are expected to break this degeneracy at non-linear scales 
even when P(k)b 2 is fixed: from Eq. 2.21-2.22, we see that / and b do not appear always in the 
same combination in the A and B functions. In particular A can be expressed as b 3 A(k,[i,/3) and 
B as b A B{k fi, j3). Since in this paper the input is the dark matter power spectrum, the degeneracy 
between b and / is not perfect, and there is the possibility of recovering these parameters separately 
with certain accuracy even for the Kaiser and Scoccimarro models. For simplicity, we do not model 
or fit the scale dependence of the bias function b(k), but instead assume that we do know the scale 
dependence of the bias, and try to recover the growth rate / and the bias amplitude, At, defined to 
be 1, 

b(k) -> A b b(k). (4.19) 

Realistic approaches need an analytical modeling of the scale dependence of the bias, that in principle 
one could expand perturbatively e.g., [52, 67]. Here we focus on the modeling of the dark matter and 
redshift space distortions and thus we measure the bias directly from the simulations. 

First, we assume that we also know that At, = 1 when we fit /. The fitting results are shown in 
Fig. 9 using the same color and line notation as in Fig. 8. As before, the error bars correspond to 
the interval defined by Ax 2 =1.0. 

As we have commented for the dark matter case, these errors are much smaller than the ones 
which could be obtained form a real survey. They contain information of 160 realizations of a large 
volume, the total volume is close to that enclosed by an all-sky survey up to z = 25. 

For both z = and z = 0.5 the models that do best at recovering / are TNS + 2L-SPT (solid blue 
line) or TNS + 2L-RPT-A2 (solid orange line), especially the latter one. In particular, at Z = 0, the 
TNS model with 2L-RPT-A<2 implementation is able to recover the / parameter within a 5% accuracy 
up to scales k = 0.17ft./Mpc and up to scales k = 0.13/i/Mpc at z = 0.5. Scoccimarro models with 
accurate implementations of the real-space power spectra such as 2L-SPT and 2L-RPT-A2 largely 
over-predict / at k > 0.05 h/Mpc. 

Next we show in Fig. 10 the results of extracting both At, and / from the halo monopole at the 
same time. We only show the two models that performed best at extracting / when At, is known: 
TNS + 2L-SPT (solid blue line) and TNS + 2L-RPT-Af 2 (solid orange line). Since now there are 
two free parameters, the 1 — a error bars correspond to contours of A% 2 = 2.3 in the f — b space. 
We see that in general, at z = and z = 0.5, 2L-RPT-A/2 works better than 2L-SPT, as we have 
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Figure 9. Best-fit values for / derived from the halo monopole for z = (left panel) and z — 0.5 (right 
panel), assuming that the bias amplitude Ab is perfectly known. Results from different theoretical models are 
shown: dashed lines correspond to Scoccimarro model and solid lines to the TNS model. Different PT models 
are shown: 1L-SPT (red lines), 2L-SPT (blue lines), 2L-RPT-M (green lines) and 2L-RPT-A/" 2 (orange lines) 
whereas 5% deviations are shown in the horizontal black dot-dashed lines, as indicated. Errors correspond to 
1 — a error, or A% 2 = 1.0. 




0.15 

. max [h/Mpc] 




Figure 10. Estimates for the parameter / corresponding to halo monopole for z = (left panel) and z = 0.5 
(right panel), treating / and Ab as free parameters. Only the models TNS + 2L-SPT (solid blue line) and 
TNS + 2L-RPT-A/2 (solid orange line) are shown for simplicity. Errors correspond to 1 — a, or A\ 2 = 2.3. 



already observed in Fig. 9. However, 2L-RPT-.A/2 tends to overestimate the bias amplitude and 
underestimate the / parameter. In particular, at z = and z = 0.5, at scales around k ~ 0.11 h/Mpc, 
this model overestimate A b by 1% and underestimate / by about 15%. However, when Aj, was set 
1, the systematic error of / was less than 5%. Similar results were obtained by [12], when the TNS 
model was combined with the Closure perturbation theory of [53]. We have checked that this large 
underestimation in / is due to the 1% overestimation in the bias amplitude. In other words, if the bias 
amplitude is set by hand to a fixed value of A b — 1.01, we obtain similar results as in Fig 10. Thus, 
this formalism tends to underestimate / by ~ 15% while overestimating A b by only < 1%. Since the 
Scoccimarro model + 2-loop perturbation theory predictions were not able to predict / when At, was 
assumed to be 1, the TNS terms A and B are the key ingredient of the TNS model to achieve a high 
accuracy recovering / when the bias amplitude is assumed. However, the signal-to-noise ratio of these 
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Figure 11. Top subpanels: different contribution terms to the halo-monopole for 2L-RPT-A/2 model: Scocci- 
marro model (red line), TNS model (blue model), TNS-A function (green line), TNS-B function (orange line) 
and P ee (black line). Bottom subpanels: ratio between these models and P tt , or the signal-to-noise (S/N). 
Dashed lines indicate a negative contribution. Horizontal dotted lines in bottom subpanels mark the reference 
quantities: S/N=l, 0.5, and 0.1. Left panels at z = and right panels at z = 0.5. 

terms is low compared to the signal-to- noise of the whole monopole term. In Fig. 11 we show this for 
z = (left panel) and z = 0.5 (right panel). In top subpanels, the total contribution of Scoccimarro 
and TNS models is shown in red and blue lines respectively. In green and orange lines, the isolated 
contribution of A and B of the TNS model is shown respectively. Black line shows the stochastic 
noise, P ee . In bottom subpanels, the ratio of all the signal terms with respect to P cc is shown. As a 
reference, the lines where S/N = 1, 0.5 and 0.1 are also shown as black dotted lines. We see that the 
signal associated to the A and B terms is much less than for the Pgs, P&e and Pgg term. In particular, 
we see that for the scales of interest (0.10 h/Mpc < k < 0.15 h/Mpc) at z — the signal for A and B 
terms is about S/N = 0.1 and even less at z = 0.5. This means that in Fig. 9, the crucial difference 
between red and blue lines (for a given PT model) comes from terms with low signal-to-noise ratio. 
In Fig. 10, when / and At, are allowed to vary, the low signal that the A and B terms have, has to 
be split to find one more parameter, and then, the accuracy recovering / must decrease. In order 
to break the degeneracy between / and b from redshift space distortions, the signal-to-noise of the 
non-linear A and B terms must be as high as possible and it could be optimized by selecting haloes 
with suitable cuts in mass. The shot noise increases with the mass cut but also does the bias. When 
optimizing for a measurement of the angle-averaged power spectrum, for a fixed P mm , the signal-to- 
noise scales as a function of mass cut like 6 2 /shot noise. If the shot noise can be approximated as 
Poissonian then the signal-to-noise can be roughly approximated as b 2 n with n number density of 
tracers. When considering redshift space distortions, the signal for the A and B terms scales like b 3 
and b A suggesting a different scaling of the signal-to-noise with mass than for the real-space power 
spectrum, which is favored by higher bias. These considerations might be useful when optimizing a 
survey selection of targets, although the bias of observable tracers might not behave as the halo bias 
especially when multiple galaxies occupy the same halo. 

5 Summary and conclusions 

Using a suite of 160 N-body simulations each with a volume of Vt, x — 13.8 (Gpc/h) 3 , we have inves- 
tigated the accuracy of analytic models in predicting the non-linear power spectrum of matter and 
dark-matter haloes in real and redshift space. The total simulated volume amounts to 2,200 (Gpc/h) 3 , 
much larger than the volume surveyed by any forthcoming or planned survey ensuring that statistical 
errors in the determination of the simulation data points is negligible. We make the dark matter and 
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halo power spectra publicly available and also provide the multipoles measured from these simulations 
for possible comparisons 7 . 

We considered a number of theoretical schemes obtained by combining standard or resummed 
perturbation theory with analytical models for redshift-space distortions (based on the simplification 
of Eq. 2.7). To predict power spectra in real space, we have employed 1- and 2-loop standard pertur- 
bation theory, and the resummed perturbation theory proposed by [24, 56] that we have generalized 
to account for 2-loop correction terms in the resummed propagator. For the redshift-space power 
spectra, we have focused on the models proposed by Kaiser [26], Scoccimarro [27] and Taruya et al. 
[28]. 

At the level of dark matter in real space, increasing the order in loop corrections for the resummed 
propagator improves the theoretical predictions of the power spectrum. In particular, working at 2- 
loop correction in the resummed propagator, A/2, an accuracy of < 1% is achieved at different redshifts 
up to the following scales: k = 0.10/i/Mpc at z = 0; k = 0.15/i/Mpc at z = 0.5; k = 0.20/i/Mpc 
at z = 1.0 and k = 0.25/i/Mpc at z = 1.5 for a 2-loop truncation of the infinite series. In general, 
working at 2-loop correction in the resummed propagator provides a more accurate description than 
working at 1-loop correction (as many of the public codes do [50, 54, 56]). 

Also, the price of working at 2-loop instead of 1-loop correction in the resummed propagator 
is not very high in terms of computational resources. It is true that evaluating A/2 involves the 
5-dimensional integration of P15, but the angular part of this function can be either analytically 
computed, or numerically precomputed for any shape of the linear power spectrum so, in the end, one 
ends up with a 2-dimensional integration, which can be easily performed. 

For dark matter in redshift space, our results show that the model by Taruya et al. combined 
with a Lorentzian damping term for the FoG effect with ero as a free parameter, is able to reproduce the 
multipoles from N-body simulations with high accuracy. In this paper, we have fit cto as a function of 
fc max using monopole, quadrupole and hexadecapole data separately, which yields to 3 different values 
for er . Although is possible to fit all these 3 multipoles with the same value of <To, the accuracy is 
expected to be reduced. When 3 different <7o are used, the high accuracy holds true for the monopole 
(< 2%) and the quadrupole (< 5%) irrespective of the flavor of perturbation theory adopted to 
compute the real-space power spectrum. For higher-order moments, such as the hexadecapole, the 
level of accuracy is more modest, ~ 10% at z = 1. This suggests that even on mildly non-linear scales 
at redshifts z < 1.5 the accuracy of the modeling of the redshift space distortions is more important 
than the modeling of the non-linear evolution of the real-space dark matter power spectrum. 

The difference between the Taruya et al. model (which attempts to include the density- velocity 
coupling at higher orders) and the other two models become more evident as the multipole order is 
increased, possibly suggesting that non-linearities become more important for higher-order multipoles. 

While in linear theory only the monopole, quadruple and hexadecapole are non-zero, in principle 
non-linearities should "excite" all higher-order multipoles, and thus cosmological signal could, in 
principle, be extracted from them. We find from the N-body simulations that the signal-to-noise 
decreases with increasing multipole, making the hexadecapole errors large even from the large suite 
of simulations considered here. This suggests that for cosmological applications most of the signal- 
to-noise at these mildly non-linear scales -where analytic approaches can provide a good modeling- 
is still enclosed in the monopole and quadrupole. 

The imprint of the baryon acoustic oscillations, is also visible in the multipoles. We find that 
all models of redshift space distortions considered do not bias the wiggles location although the more 
linear models (Kaiser and Scoccimarro) over-predict their amplitude. These considerations might be 
relevant for recovering in an unbiased way the angular and radial BAO information (separately) from 
forthcoming surveys. 

Overall, the accuracy of the analytic description allows measurements of the logarithmic growth 
rate / to percent level, when 00 is allowed to vary. In this case, when the Taruya et al. model 
is combined with the RPT prediction for the power spectrum in real space, / is recovered within 
< 1% up to fc max = 0.15/i/Mpc for the monopole and quadrupole (separately) at z = and up to 

7 http : //ice .ub . edu/~hector/Hector_Gil_Marin/Public .html 
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k max = 0.20ft-/Mpc at z = 1. This indicates that the Taruya et al. model combined with the RPT 
real-space predictions is accurate enough to be used for precision cosmology. 

Since most of the current and future redshift surveys target galaxies as tracers of the matter 
distribution, a more realistic way of estimating / is to use dark-matter haloes instead of the dark 
matter density. The limited mass resolution of our N-body simulations, allows us to consider only 
cluster-sized haloes M > 1O 14 M / h. Dealing with isolated haloes has the advantage of eliminating the 
imprint of FoGs from the power spectrum. However, the effect of the scale-dependent (and possibly 
non-linear) bias plays an important role. In this work, we have assumed that the bias is linear and 
that its /c-dependence is known. Under these approximations, we have been able to recover / with 
< 5% when the amplitude of the bias is known a priori. For these massive haloes, the effect of bias 
is important: the degraded accuracy in recovering / indicates that, at least for these massive haloes, 
the modeling of biasing is crucial. In particular, given the high shot noise that the statistics of these 
tracers have, a modeling of its behavior both as a function of scale and halo mass is important at 
mildly non- linear scales. 

When both the / parameter and the bias amplitude are allowed to vary, we recover the bias 
amplitude to 1% level in the best cases, but / is underestimated by 10% — 20% at z — and slightly 
more at z — 0.5. Similar results are reported in the literature [12]. Most likely this is not (or not only) 
due to a limitation of the model for the power spectrum, but also to the poor signal-to-noise ratio 
of the population of haloes used to extract /. It remains to be seen whether reducing the halo mass 
threshold increases the signal-to-noise. The scaling with halo mass of the signal-to-noise of the RSD 
terms that can break the degeneracy between / and b, is different from the usual nP ~ b 2 n used for 
the angle-averaged, real-space power spectrum. Our simple considerations indicate that a lower mass 
threshold increases the signal-to-noise, but if the number density is kept fixed, then the signal-to-noise 
for the A and B terms is favored by a higher bias (the signal-to-noise scales like b 3 n and b 4 n rather 
than like b 2 n). These considerations might be useful when optimizing a survey selection of targets, 
although one should keep in mind that bias of observable tracers might not behave as the bias of the 
host haloes. 

The results found in this paper are in agreement with those found in recent works. Kwan et al. 
[12] found that, most of the RSD models fail at recovering /, underestimating its value even at large 
scales for z = and z = 0.5. Considering the Taruya et al. model, relaxes the discrepancy (compared 
to Kaiser or Scoccimarro models), but does not completely fixes the problem. Other works, such as 
de la Torre & Guzzo [11], show that the Taruya et al. model (in configuration space) + numerical 
and phenomenological schemes to estimate the real-space spectra, are able to recover / with < 5% 
accuracy from different low-biased galaxy population. Okumura & Jing [14] show that (3 can be 
estimated accurately from very massive haloes {Mh > 1Q 1a Mq/K) using linear theory (Kaiser model) 
at large scales (k < 0.1/i/Mpc); and / can be also recovered with similar accuracy when the bias is 
assumed to be known, although the precision achieved is not very high. 

Future surveys and missions will provide datasets about the distribution of galaxies on large 
scales. We envision that in order to extract useful information from these datasets, we will need more 
accurate theoretical models of structure formation. In this paper, we have shown that for extracting 
the growth of structure correctly, accurate analytic models for both real and redshift space clustering 
are crucial. In particular 2L-RPT-A2 in combination with the Taruya et al. formula, seems to be able 
to recover / accurately from massive haloes, when the bias is assumed to be known or equivalently 
when one wants to recover the combination j3 — / jb. However, in our study we found that this 
approach fails when trying to recover the individual values of / and b simultaneously. This can be 
due to the low signal-to-noise ratio of this halo population, or a limitation in the model itself. Should 
the determination of / and b separately rather than in the /3 combination from RSD alone become a 
priority, studies will attempt to improve both the real and redshift-space models. In the case of real 
space, this can be done extending resummation theories both in Eulerian or Lagrangian spaces. In the 
case of redshift space, considering higher-order terms in the TNS formula should improve the model. 
Finally, it is also important to model correctly the stochasticity associated with the halo population 
and to determine correctly the scale dependence of the bias and its possible non-linearities. 
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A Standard Perturbation Theory 

In this section, we provide a short summary of the equations used to compute the 1- and 2-loop 
correction in Eulerian perturbation theory. For a detailed description of perturbation theory see 
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[31, 68]. 

According to standard perturbation theory (SPT) the power spectrum in real space can be 
expressed as a sum of loop corrections, 

P(k) = p(°) (Jfe) + P« (k) + P (2) (*)..., (A.l) 

where P^°'(k) = P lln (fc) is the linear term. For Gaussian initial conditions, the different loop correc- 
tions read as, 

P (1) (/c) = 2P 13 (k) + P 22 (k) 1-loop correction, (A.2) 
P [2) {k) = 2P 15 (k) + 2P 24 (k) + P 33 (k) 2-loop correction, (A.3) 

where, as already mentioned in the main text, the subscripts i and j refer to the perturbative order of 
the terms <5(k) used in eq. 2.1 to compute the power spectrum Py(fc). For the case of matter-matter 
power spectrum, namely Pss, these terms are [32], 



P 13 (fc) = 3P lin (fc) / ^PKk^-q)^), (A.4) 



lin (fc) J "' C ' i7 ^ ' - - ■ n '"" • 

d 3 q 

(2^)3 



(2tt) e 

Ps>{h = 2 / ,^ :1 F| 2 (q,k-q)/> K M ). (A.r,l 
d 3 qi d 3 q 2 T , _ ,„,„,, . ,> •. ^6) 



P 15 (fc) = 15P lin (fc) J ^3 ( 2^3 f 5 qi, -qi, qa, -q 2 )P lin (gi)P lin ( g2 ) 

^ 24 (fc) = 12 y ^^P 2 s (q 1 ,k-q 1 )P|(q 1 ,k-q 1 , cu, -q 2 ) x (A.7) 
x P lin (q 1 )P lin ( (Z2 )P lin (|k-q 1 |), 

[/ ^P 3 s (k,q,-q)P lin (<z) n 



P 33 (fc) = 9P lin (fc) 



(A. 



6 J (2^3 ^ ^I 2 (qi ; q2,k- qi -q 2 )P lin (gi)P lin (g 2 )P lin (|k-q 1 -q 2 [). 



In the 1-loop correction, P 22 accounts for the mode coupling between vectors with frequencies k q 
and q, whereas Pi 3 can be interpreted as the 1-loop correction to the linear propagator. In a similar 
way, in the 2-loop correction term, only the second term of P33 accounts for a full 2-loop mode coupling 
because is the only term that contains P lm (|k — qi — q 2 |). Also note that P 2 4 contains a term similar 
to a 1-loop coupling, P lm (|k — qi |), is similar to P 22 . Pi 5 and the first term of P33 contain no coupling 
between k and q.; and can be interpreted as a 2-loop propagators. In particular, the full n-propagator 
can be written as [32] 

Pi n (fe) = n!!P lin (fc) / 5% • • • §^ F n( k > Qi. "Qi. • • • . q„ -q.)^ lin (9i) • ■ ■ ^"(fe), (A.9) 
where x = (n — l)/2. 

These similarities between these terms is the basis of the resummation process that is described 
in Appendix B. The kernels of Eqs. A.4-A.8 are expressed as, 

P n (qi,...,q n ) = V ^ m( ^;; ; • ' q? " } [(2n + l)a(k, k 1 )F n _ m (q m+1 , . . . , q n )+ (A.10) 

+ 2/3(k,ki,k 2 )G„_ m (q m+ i, . . . ,q„)] , 
n G ( 1 

G n (qi,...,q n ) = Y] m . U "\ [3a(k,ki)P»_ m (q m+ i,...,q„)+ (ATI) 
^ (2n + 3)(ra- 1) 

+ 2n/3(k,ki,k 2 )G„_ m (q m+ i, . . . ,q n )] , 
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with Fi = 1 and G\ = 1. Also, ki = qi H + q m , k 2 = q m +i + • • • + q„, k = ki + k 2 where the 

functions a and /3 are defined as, 

a(k,ki) = (A.12) 
/3(k,k 1; k 2 )^ fc2 j, k / 2 , 2 k2) . (A.13) 
The symmetrization process of the kernels is given by, 

■fn(qi) ■ ■ • ,q„) = -f ^^"(^(l)' • • • '^(n))' ( A - 14 ) 

7T 

where the sum is taken over all the permutations 7r of the set {1, . . . , n}. In particular, the expressions 
for F|(ki,k 2 ) and G 2 (ki,k 2 ) are, 

F| (k x ,k 2 ) = 7 + ^os^| + + ?cos 2 0, (A.15) 

G^k l5 k 2 ) = | + icos0(| + ~)+~cos 2 0, (A.16) 

where, cos = (ki ■ k 2 )/(fciA; 2 ). 

This SPT formalism presents some shortcomings. As noted by [24], at large scales only the linear 
term contributes to the total power spectrum. However, at smaller scales, all loop corrections become 
of the same order with a significant cancellation among them. In particular, at low redshifts it can be 
seen that 1-loop correction overestimates the full power spectrum (from N-body simulations), whereas 
the 2-loop correction underestimates it. This is due to the fact that P^ is negative on large scales 
while P^ 2 ) is positive, and both almost cancel out giving a remaining quantity which is close to the full 
power spectrum. In the same way, as we go to higher order, more cancellations come out among the 
different loop corrections. Thus, truncating at certain loop in SPT will naturally produce a systematic 
over- and under-prediction of the real-space power spectrum. A way to avoid this behavior is to resum 
some terms of the total SPT expansion. In [24, 25] a formalism for resumming part of these terms 
was proposed. The resulting expansion presents a more controlled behavior because each different 
loop contributes only positively to the total power spectrum and acts at different scales. In Appendix 
B we present an alternative description (but mathematically identical) of the resummation presented 
by [24, 25]. As an extension of current works, we write not only the 1-loop resummed propagator, 
but we perform our computation up to 3-loops. 



B Resummation in Standard Perturbation Theory 

In this section we present the derivation of Eq. 2.11, 2.14 and 2.15. We also show that the Afi 
expression is identical to the one used in [56] when the propagator is perturbed at 1-loop. For 
completeness, we also show the result of perturbing the propagator for 3-loops. These equations come 
from resumming some terms in SPT under certain approximation in the kernels. In particular, under 
the Zel'dovich approximation, the kernels read [69], 



where k = ki 



■ ) k n ) 



G£(ki, ...,k„) 



1 k • ki k • k„ 



1 k • ki k • k n 



(B.l) 
(B.2) 



k„ . As shown by [24] , with this approximation the resummation process yields, 



P(k, z) - P lin (k, z) + P£(k, z) + P 3 2 3 L (fc, z) 



p. 



(n-l)L 



(k,z) 



M Q {k,zf 



(B.3) 
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where, 



with 



JVo(k, z) = exp 



1 



k 2 a v {zf 



(B.4) 



n-l 



x P 11 



(2tt) 3 (2tt) 2 



Z) X 



(B.5) 



and cr„ is a characteristic scale defined as 



4tt 



dq 



P lin (q,z) 



(B.6) 



We will refer to Eq. B.3 as RPT-Ao model. With this technique the behavior of PT improves, since 
every new loop adds a positive term that only acts on a small range of scales. Therefore, using this 
technique the oscillatory behavior observed in standard PT vanishes. 

In this section we show that performing a slightly different approximation (not Zel'dovich) in 
the kernels we can end up with the same formula used in [56]. Also, depending on how we 'factorizc' 
the kernels, we will end up with a 1-, 2- or higher-loop correction in the resummed propagator. These 
formulae present a notable improvement respect to Eq. B.3 for a similar computational effort. In 
this work we do not follow the approach of Feynman diagrams to resum the infinite terms as it is 
done in [24, 25]. Alternatively, we present a different approach for doing this without requiring any 
knowledge of quantum field theory. We hope that this way of resumming is clearer for the reader who 
is no familiar with this kind of formalism. Furthermore, this approach allows us to easily compute the 
resummed propagator for higher-order loops. Our method consists of rewriting the terms of the £-loop 
correction (where I = (n + m)/2 — 1), namely P nm , as a sum of subterms which can be associated to 
lower loop corrections as we show below 8 , 



Pnm(k) — Pnm(k) + Pnm(^) + 



2L 

nm 



(k) + ■ ■ ■ + P^ n {k) 



(B.7) 



The subterm with an index OL contains a P lln (fc) and corresponds to the linear power spectrum (no- 
loop correction), the subterm with an index 1L contains a P lln (|k — qi|) and therefore is similar to 
p22(fc) (that corresponds to 1-loop correction). In the same way, the 2L subterm is similar to p33(fc) 
(that corresponds to 2-loop correction) because contains a term P lln (|k — qi — q 2 |) and so on. The 
generic way of writing these terms is the following 9 . 
The 0-L subterm can be written, 



1. for n & m even, 



Pnm(k) — 0; 



(B.8) 



2. for n & to odd, 



P^(k) = n!!m!! 



(2tt) 3 ' ' 

x P£(k, q r,-qr,- 



(2.)3 (27r) 3 - (27r) 3^(k,q?,-q?,..,<,-qSJ X (B.9) 
. , q™ , -q™ )P lin (fc)P lin (^) ■ ■ • P Un (C )P lin (9i m ) • ■ • P Un (qZ)> 



where Xi = (i — l)/2. 
The 1-L subterm can be written, 



8 The redshift dependence is understood for simplicity: it only appears through P lln . 

9 Note that the terms with n odd and m even (and vice versa) vanish for Gaussian initial conditions. 
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1. for n & m even, 

P%kW = £n(»-l)!!m(m-l)!! 



d 3 qi d 3 q% d 3 q™ n d 3 q? 



(2tt) 3 (2tt) 3 ' ' ' (2tt) 3 (2vr) 3 ' ' ' (2tt) 3 



j3 m 

'" m P lin (|k- qi |)P lin ( gi ) x 



x F»( qi) k - qi! q", -q* . . . , q" , -q« )P lin (g 2 ") . . . P lln (<z" ) x 



ilin / 



x P4(qi,k-qi,q™ 
where Xi = z/2; 
2. for n & to odd, 

The 2-L subterms can be written, 

1. for n & to even, 

2. for n & to odd, 
n!!™!! 



-q 2 



q™ )P lm (gD...P lin (g™ 



(B.10) 



p2L _ q. 



(B.ll) 



(B.12) 



P 



2L 



6 



(n- l)(m- 1) 



d 3 9l d 3 g 2 d 3 gj d 3 q n x d 3 ql n d 3 q^ 



(2tt) 3 (2tt) 3 (2tt) 3 ' ' ' (2tt) 3 (2tt) 3 ' ' ' (2tt) 3 



P lm ( qi )P lm (q 2 ) x 



X P*( qi , q 2 , k Ql q 2 , q£, -eg, . . . , q^ , -q^)P lln (|k - q x - q 2 |)P lm (<Z3 n ) ■ ■ ■ i^feSJ 



x -F , m(qi,q2,k-q 1 -q2,q£\-qS\...,q£ 

with = (i — 3)/2. 
The 3-L subterms are, 
1. for n & to even, 

^(*) = ^"(« - 2)(n - l)!!m(m - 2)(m - 1)!! 



q ™ )P 1 ^ 3 m )---^ lm (C l ), 



lin/ 



(B.13) 



d 3 gi d 3 g 2 d 3 g 3 A" rf 3 C 



(2tt) 3 (2tt) 3 (2tt) 3 (2tt) 3 ' ' ' (2tt) 3 



d 3 qT d 3 q™ 



^(qi,q2,q3,k- qi - q 2 - q 3 ,q4,-q?: 



-q^-qij x 



(2tt) 3 ■ ' ■ (2tt) 3 

x Pm(qi>q2,q3,k-qi -q 2 -q3,q™,-q™,---,qr„ 

x P lhl (g 3 )P lin (|k - qi - q 2 - q 3 |)P lin (gJ) ■ • • P^{^L)P^(.Qa) ■ ■ ■ ^(CJ 



with Xi = (i - 4) /2; 
2. for n & to odd, 



q™ )p lin (gi)p lm (g 2 ) X 

(B.14) 
(B.15) 



and similarly for higher-order subterms. In particular, it is important to note that when n = m the 
subterm with the highest loop correction (I = n — 1) is expressed as, 



Pt- 1)h (k) = n\ 



d 3 qi d 3 q2 d 3 q n -\ 2 \^ \ D lin/ \ Dlin/ \ 

-^r F n (qiiqa,--- ,qn-i,k- 2_^m)P l m 0i).-. p 



x P 



lin/ 



(2tt) 3 (2tt) 3 '■■ (2tt) 3 

n-1 



(B.16) 



which is the same term used in Eq. B.3. This indicates that the resummation of terms described in 
[24] corresponds to resumming the terms P„„ for I < (n + to)/2 — 2. In order to make possible the 
resummation we perform an approximation in the kernels. 
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B.l 1-loop factorization 

If we want to end up with a resummed propagator of 1-loop correction, the prescription in factorizing 
the kernels is the following, 

1. for 0-L subtcrms, 



F' (k, q X) -qi, . . . , q^ , -cb B ) ~ i [3!F 3 S (k, qi, -qi)] • . . . ■ [3!F 3 S (k, q Xn , -q Xn )] ; (B.17) 

2. for 1-L subterms, 

^(qi.k-qi.qa.-qa,...,^,-^) ~ ~ [2LF| (q X) k - q x )] x (B.18) 

x [3LF 3 s (k,q 2 ,-q 2 )] ^'(k^.-cbj]; 

3. for 2-L subterms, 

^(qi,q 2 ,k- qi - q 2 ,q 3 ,-q 3 , ■ ■ ■,<!*„,, -QxJ - ~, [3!P 3 s (qi, q 2 , k - q x - q 2 )] X (B.19) 

n! 

x [3!F|(k, q 3 , -q 3 )] • . . . • [3!F|(k, q^, -q x J] ; 

and similarly for higher-order loops. We factorize the n-order kernel in a product of 2- and 3-order 
kernels (which correspond to 1-loop correction terms) keeping the sum qi + q 2 + • • • = k in all the 2- 
and 3-order kernels. Under these approximations we rewrite the subterms of Eq. B.7. 

1. For 0-L subterms with n & m odd, 

P ™ {k) -x n \x m l\P**(k)) ' (B ' 20) 

with Xi = (i — l)/2. 

2. For 1-L subterms with n & m even, 

with Xi = {i- 2)/2. 

3. For 2-L subterms with n & m odd, 

P nrn(k) 3^ ( ) (B.22) 

with Xi = (i — 3)/2; 



x n \x m \ \P l ™(k) 



and similarly for higher-order subterms. Now we can proceed with the resummation of Eq. B.7. 
Reordering the terms we write, 



odd odd even even odd odd 



P W = E E P nt(k) + E E P nrn(k) + E E P ™( k ) +.... (B.23) 



n— 1 m—1 n=2 m— 2 n— 3 m— 3 

The first term is, 



odd odd oo oo 

EE«=^)E E ^ 

n—1 m—1 x n —0 Xm—O 



P Vm (k) 



P lin (fc)exp[2P 13 (fc)/P lin (fc)] . 

(B.24) 
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The second term is, 



even even 



oo oo 



E E p n h m (k) = p 22 (k) e E 



The third term is, 

odd odd 



x n =Q x m =0 



OO OO 



EE^w-^wE E 



1 



.!.*_! 



n— 3 m — 3 



x n =Q x m =0 



Piajk) 
P lin (k) 



Pi 3 (k) 



X-n+Xn 



= P 22 (fc)exp[2P 13 (fc)/P lin (fc)] 



(B.25) 



X n -\-X n 



= P 3 2 3 L (fc)exp[2P 13 (£0/P lin (fc)] 



_P iin {k)_ 

(B.26) 

and the same for higher-order terms. Therefore, after this resummation we can express the power 
spectrum as, 

P(k) = \P nn (k) + P£(k) + P 2 3 L (fc) + ...} M(*) a , (B.27) 



with, 



M(fc) = e 3 q>[^3(fc)/P Un (*0] 



(B.28) 



Here, P^CO = P22(fc) is given by Eq A. 5, whereas P 3 2 3 L (£;) is given by the second term of Eq. A. 8, 



d 3 gi d 3 g 2 



) 3 (2tt) ; 



Pf (q 1; q 2 , k - qi - q 2 )P lm ( qi )P lm (q 2 )P^(\k - qi - q 2 



(B.29) 



The factor in the exponential, Pi 3 (k)/P(k), can be partially computed, because the integral over the 
angular part of Eq. A. 4 can be performed for any shape of P lm . For doing this, is convenient to 
express the symmetrized P| kernel as, 



6F3 S (k,q,-q) = ^(k, q) [7a(k, k + q) + 4/3(k, k + q, -q) 



(B.30) 



G s 2 (k, -q) [7«(k, k - q) + 4/?(k, k - q, q) 



+ - a (k,q)[P|(k,-q)-P|(k,q)]. 

Taking into account Eqs. A.12-A.15, the angular dependence is now explicit through cos 8, and the 
angular integration of Eq. A. 4 can be performed. This computation yields, 



, Pisjk) = r 4tt 
'P lin (fc) i 504fc 3 q 



- 79fc 5 g 3 + 50g 5 fc 3 - 21kq 7 - 



(B.31) 



+ ^(fc 2 -<z 2 ) 3 (2fc 2 + 7g 2 )ln 



k — q 



k + q 



P lin (q)dq, 



which is the expression proposed by [56] as the resummed propagator. We will refer to the model of 
Eq. B.27 as RPT-Wi model. 

B.2 2-loop factorization 

It is also possible to obtain a 2-loop resummed propagator if we factorize the kernels in a different 
way. In that case, we want to split the kernels in pieces of P 3 S and P| kernels, 

1. for 0-L subterms with n = 1, 5, 9, 13, . . . , 



^n(k,qi,-qi, ... ,qx„,-qa„) - [5!P 5 s (k, qi, -qi, q 2 , -q 2 )] 
x ... • [5!F 5 *(k,q a;n _i,-q a!n _i ) q a!n ,-q a!n )] ; 



(B.32) 
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2. for 0-L subterms with n = 3, 7, 11, 15, ... , 



^(k,qi,-qi, 



1 



1 HX n , 4l„ 



[3!P|(k, qi , -qj) • 5!P 5 s (k, q 2 , -q,, q 3 , -q 3 )] (B.33) 



x [&Fg(k,ci Xn -i,-q Xn -i,q Xn ,-q Xn )}; 
3. for 1-L subterms with n — 4,8, 12, 16, ... , 

F n s ( qi ,k-q l!q2 ,-q 2 , ... ,q ;c , i ,-q :c J~i[2!F 2 5 (q 1 ,k-q 1 )].[3!F3 S (k,q 2 ,-q 2 )]x (B.34) 



x [5!i ? g( k >q3,-q3,q4,-q4)] 

4. for 1-L subterms with n = 2, 6, 10, 14, ... , 

-f^(qi,k- qi,q 2 ,-q 2 , ... ,q x „,-q x J~ — [2!P|(qi,k - q x )] x 



[5!F 5 S (k, q x „_i, -qx n -i, q s „, -q^ 



n! 

x [5!P 5 s (k,q 2 ,-q 2 ,q 3 ,-q 3 )] 
5. for 2-L subterms with n = 1, 5, 9, 13, . . . , 

1 



(B.35) 



[5!P|(k, q a „_i, -q^-i, q Xn , -q^J] 



^(qi,q2,k-qi -q 2 ,q 3 ,-q3, •■• ,q*„,-qz„) - -r [3!F 3 s (qi,q 2! k- qi - q 2 )] • [3!P 3 s (k, q 3 , -q 3 )] x 



x [5!F 5 s (k,q 4 ,-q4,q5,-q5)] • • • [5!F|(k, q x „_i, -q K „-i, q^ , -QxJ] ; 

(B.36) 

6. for 2-L subterms with n = 3,7, 11, 15, ... , 

i^(qi,q2,k- qi - q 2 ,q 3 ,-q 3 , ... , q Xn , -q x J ~ ^ [3!P 3 (qi, q 2 , k - q x - q 2 )] x (B.37) 

x [5!F 5 s (k,q 3 ,-q 3 ,q 4 ,-q 4 )] • ... ■ [5!F 5 s (k, q x „_i, -q^-i, q Xn , -q*J] ; 

and similar for higher-order loops. As before, with this approximation we can perform an exact 
resummation of all terms yielding to, 



P(k) = [P bn (fc) + P%{k) + P$ 3 L (k) + ...} N 2 {k) 



(B.38) 



with, 



N 2 {k) = cosh 



>2P 15 (k) 
P lin (fc) 



+ 



Piz(k) P' ! "(fc) . 
P»»(fc) V 2P 15 (fc) 



sinh 



'2P 15 (k) 
P lin (fc) 



(B.39) 



This is the general form for J\f 2 - However, at large scales, Pi 5 < 0, and this expression becomes, 
J\f 2 (k) = cos 



'2|Pis(fc)| 



P lin (fc) 



Pi 3 (fc) / P n "(fc) . 
P li «(fc)V2|Pi5(fc)| Sln 



/2|Pi B (*)| 



P lin (fc) 



(B.40) 



We will refer to Eq. B.38 as RPT-A/2 model. 

Note that if we perform on Pi5(fc) in Eq. B.39 the approximation, 



P|(k, qi, -qi, q 2 , -q 2 ) ~ ^3!P|(k, qi , -qi)3!F|(k, q 2 , -q 2 ), 



(B.41) 



we obtain that 2Pi 5 — > (P i3 /P lln ) 2 P lln and therefore J\f 2 — > jVi . In the same way, when we apply the 
Zel'dovich approximation on the kernel of Pi 3 , 



P 3 s (k,q,-q) 



1 k • q — k • q 

3!~^ ' 



(B.42) 
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we obtain that 2Pi3/P lm — > — k 2 a 2 , and therefore, A/i — > No. Note that for the computation of Ni(k) 
one needs to compute Pis(fc) which requires the knowledge of P|. This computation is a 6-dimensional 
integral that reduces trivially to 5-dimensional exploiting rotational invariance. In principle, one could 
integrate analytically the remaining 3 angles and reduce the computation of Pi 5 to a 2-dimensional 
integral in the same way Pi 3 is reduced to a 1-dimensional integral in Eq. B.31. However this is 
hard, because the symmetrized kernel Fg is the sum of 5! = 120 different cyclic permutations. A 
possible alternative, is to precompute the angular part of Pi 5 as a 3-dimensional integral for a wide 
range values of k, q\ and q%, and then use this to compute Pi 5 as a 2-dimensional integral for any 
shape of P lm . Nevertheless, for practical reasons, in this paper P i5 is always computed numerically 
as a 5-dimensional integral. For completeness we also report the expression for the 3-loop resummcd 
propagator, A3 (k) function as function of the full-propagator terms P13, Pig and P17, 



N 3 (k) 



/6Pi 7 (fc) 
P lin (fc) 



P 13 (fc) 3 P*"(k) 
P lin (fc) V 6Pi 7 (fc) 



D 



lePnjk) 
P lin (fc) 



(B.43) 



2Pi 5 (fc) 3/ 
P lin (/c) V 



' P Iin (fc) 
6Pn(k) 



C 



3 / 6Pi7(fc) 

V P lin (fc) 



where the functions A,B and C are given by, 



A{x) = exp(a;) + 2exp ^— ^ cos ^"^^ 
B(x) ee exp(a;) - exp 
C[x) ee exp(a;) - exp 




V3 sin ( ~^~ x 



cos j —^~x I — V3sm I ~^~ x 



(B.44) 
(B.45) 
(B.46) 



We do not use this function in this paper, because it requires the computation of Pi 7 which is a 
8-dimcnsional integral (after exploiting rotational invariance), which goes beyond the scope of this 
paper. We leave the analysis of this function for a future work. 
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